% version for the calculation of shocks for all markets in turn
nm=8;
muu=load('mu_all3.txt');
eps=load('eps_all3.txt');
c_e=cov(eps); % covariance matrix of eps
sd=sqrt(diag(c_e));
k1=200; % total length
k=100; % for display purposes
hilos=load('hilos.txt');
bigT=length(eps(:,1));
eps_s=sort(eps);
giorni1=load('date_7_18_95.txt');
MATLABDate = x2mdate(giorni1);
giorni=datestr(MATLABDate,2);
% load matrices
aplus=load('aplus.txt');   % impact matrix lag 1 overall period dpos
aminus=load('aminus.txt'); % impact matrix lag 1 overall period dneg
acplus=load('acplus.txt');   % impact matrix lag 1 crisis period dpos
acminus=load('acminus.txt'); % impact matrix lag 1 crisis period dneg
ac2=load('a_lag2.txt'); % impact matrix lag 2 
omega1=load('om.txt');
dneg=load('dneg.txt');
beta=load('beta.txt');
ac1=beta+(aplus+aminus)/2;
ad1=(acplus+acminus)/2;
area=zeros(nm,nm);
area0=zeros(nm,nm,bigT);
area1=zeros(nm,nm,bigT);
area2=zeros(nm,nm,bigT);
spillPRE=zeros(nm,nm);
spillDURING=zeros(nm,nm);
spillPOST=zeros(nm,nm);
spillWHOLE=zeros(nm,nm);
% build the shocked solutions
for t=2:bigT
    %shk=[550 601 832 765 602 603 1274 474];
    %shk=ones(8)*1274;
    %t=shk(irf);
    % select a row from muu
    if ge(t,473) & le(t,842) % crisis period
        ac=ac1+ad1;
        om=omega1(:,1)+omega1(:,2);
    elseif ge(t,843) % post crisis period
        ac=ac1;
        om=omega1(:,1)+omega1(:,3);
    else
        ac=ac1;
        om=omega1(:,1);
    end
    mu=muu(t,:);
    mu=mu';
    nm=length(mu);
    mu=[mu ; hilos(t-1,:)'];
    a=[ac ac2; eye(nm) zeros(nm,nm)];
    om=[om ; zeros(nm,1)];
    b=zeros(nm*2,k1,nm); % difference as in IRF graphs
    b1=zeros(nm*2,k1,nm); % baseline
    b2=zeros(nm*2,k1,nm); % shocked solution
    b3=zeros(nm*2,k1,nm); % ratio
    s_mat=zeros(nm,nm);
    for mkt=1:nm
        if mkt ==1
            s=[sd(mkt) ; c_e(2:nm,mkt)/(sd(mkt))];
        elseif mkt==nm
            s=[c_e(1:nm-1,mkt)/(sd(mkt)) ; sd(mkt)];
        else
            s=[c_e(1:mkt-1,mkt)/(sd(mkt)); sd(mkt); c_e(mkt+1:nm,mkt)/(sd(mkt))];
        end
        s=zeros(nm,1);
        s(mkt)=sd(mkt);
        s_mat(:,mkt)=s;
        
        s=[s;zeros(nm,1)];
        b(:,1,mkt)=a*(s.*mu); % the vector is now multiplied by mu element by element
        b1(:,1,mkt)=om+a*mu;
        b2(:,1,mkt)=om+a*((1+s).*mu);
        for i=2:k1
            b(:,i,mkt)=a*b(:,i-1,mkt);
            b1(:,i,mkt)=om+a*b1(:,i-1,mkt);
            b2(:,i,mkt)=om+a*b2(:,i-1,mkt);
        end
    end
    b=b*log(10)*sqrt(252)*sqrt(pi)/sqrt(8);
    b1=b1*log(10)*sqrt(252)*sqrt(pi)/sqrt(8);
    b2=b2*log(10)*sqrt(252)*sqrt(pi)/sqrt(8);
    b3=b2./b1;
    area(:,:)=sum(b3(1:8,:,:),2)-k1;
    area0(:,:,t)=area;
    area1(:,:,t)=area./kron(diag(area)',ones(8,1));
    area2(:,:,t)=area./kron(sum(area,2)',ones(8,1));
end
vol_imp=mean(area0,3)
vol_imp_ij=vol_imp-diag(diag(vol_imp)); % take out the main diagonal
VSB=sum(vol_imp_ij,1)./sum(vol_imp_ij,2)'